% $Id: robpic.tex 203 2009-04-20 01:20:13Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

\magnification=1200
\input macros.tex
\input docmac.tex
\openup 1truept

\def\bmax{B_{\hbox{\ninerm MAX}}}
\def\DS{\hbox{\tt DATA\_SECTION}}
\def\PS{\hbox{\tt PARAMETER\_SECTION}}
\def\PCS{\hbox{\tt PRELIMINARY\_CALCS\_SECTION}}
\def\IS{\hbox{\tt INITIALIZATION\_SECTION}}
\def\PROS{\hbox{\tt PROCEDURE\_SECTION}}
\def\RUNS{\hbox{[RUNTIME_SECTION]}}
\def\ADM{AD Model Builder}
\def\chapno{1}
\section{Robust Nonlinear regression with \ADM}
\input macros.tex
\input docmac.tex
\input pictex
\font\rmlarge=cmr18 at 18truept
This section is intended to demonstrate the advantages 
of using \ADM's robust regression function over standard nonlinear
least square regression procedures. Further discussion about the
underlying theory can be found in the AUTODIF user's manual, but
it is not necessary to understand the theory to make use of the
procedure.
This example estimates the parameters describing a growth curve,
from a set of data consisting of ages and size-at-age data.
The form of the (von Bertalanffy) growth curve is assumed to be
$$s(a)=L_{\infty}\big(\,1-\exp(-K(a-t_0))\,\big)\eqno(.1)$$
The three parameters of the curve to be estimated are
$L_{\infty}$, $K$, and $t_0$.

Let $O_i$ and $a_i$ be the observed size and age of the $i$'th 
animal. The predicted size $s(a_i)$ is given by equation $.1$.
The least squares estimates for the parameters are found by minimizing
$$\min_{L_\infty,K,t_0} \sum_i \big(O_i -s(a_i)\,\big)^2$$.
The code for the admodel template is found in the file {\tt vonb.tpl}.
\beginexample
DATA_SECTION
  init_int nobs;
  init_matrix data(1,nobs,1,2)
  vector age(1,nobs);
  vector size(1,nobs);
PARAMETER_SECTION
  init_number Linf;
  init_number K;
  init_number t0;
  vector pred_size(1,nobs)
  objective_function_value f;
PRELIMINARY_CALCS_SECTION
  // get the data out of the columns of the data matrix 
  age=column(data,1);
  size=column(data,2);
  Linf=1.1*max(size);  // set Linf to 1.1 times the longest observed length
PROCEDURE_SECTION
  pred_size=Linf*(1.-exp(-K*(age-t0)));
  f=regression(size,pred_size);
\endexample
The only new feature introduced in this code is the use of the
{\tt regression function} which calculates the
log-likelihood function of the nonlinear regression. 
The inputs to this function are
the vector of observed quantities and the corresponding vector
of predicted quantities. 
A graph of the least-square estimated growth curve and the observed
data is given in figure~1. The parameter estimes and their estimated 
standard deviations which are produced by \ADM\ are also given.
For example the estimate for $L_\infty$ is  54.86 with a standard
deviations of 2.47. Since a 95\% confidence limit is about 
$\pm$ two standard deviations the usual 95\% confidence limit of
$L_\infty$ for this analysis would be $54.86\pm 4.94$.
\vbox{
\medskip
\quad\hbox{
\beginpicture
  \setcoordinatesystem units <.18in,.04in>
  \setplotarea x from 0 to 16.5, y from 0 to 50 
  \axis left label {Size} ticks
    numbered from 0 to 50 by 10 
  /
  \axis bottom label {Age} ticks
    numbered from 0 to 16 by 2 
  /
 \multiput {\hbox{$\bullet$}} at "vonbg.dat" 
 \put {Figure 1} at 4 45
 \plot  "vv.rep" 
\endpicture
\hfill
}
\medskip
\qquad Results for nonlinear regression with good data set
\medskip
{\openup 1pt
\beginexample
 index         value      std dev       1       2       3   
    1   Linf  5.4861e+01 2.4704e+00  1.0000
    2   K     1.7985e-01 2.7127e-02 -0.9191  1.0000
    3   t0    1.8031e-01 2.9549e-01 -0.5856  0.7821  1.0000
\endexample
}}
A disadvantage of least squares regression is the sensitivity of the 
estimates to a few ``bad'' data points or outliers. Figure~2 show the
least squares estimates when the observed size for age~2 and age~14
have been moved off the curve.
there has been a rahter large change in some of the parameter
estimates. For example the estimate for $L_\infty$ has changed
from $54.86$ to $48.91$ and the estimated standard deviation for
this parameter has increased to $5.99.$
This is a common effect of outliers on least-squares estimates.
They greatly increase  the size of the estimates of the standard deviations.
As a result the confidence limits for the parameters are increased.
In this case the 95\% confidence limits  for 
$L_\infty$ have been increased from $54.86\pm 4.94$ to $48.91\pm 11.98$.

Of course for this simple example it could be argued that a visual
examination of the residuals would identify the outliers so that
they could be removed. This is true, but in larger nonlinear models
it is often not possible or convenient to identify and remove
all the outliers in this fashion. Also the process of removing
``inconvenient'' observations from data can be uncomfortably
close to ``cooking'' the data in order to obtain the desired result
from the analysis. An alternative approach which avoids these
difficulties is to employ \ADM's robust regression procedure
which removes the undue influence of outling points without the
need to expressly remove them from the data.

\vfil
\vbox{
\medskip
\quad\hbox{
\beginpicture
  \setcoordinatesystem units <.18in,.04in>
  \setplotarea x from 0 to 16.5, y from 0 to 50 
  \axis left label {Size} ticks
    numbered from 0 to 50 by 10 
  /
  \axis bottom label {Age} ticks
    numbered from 0 to 16 by 2 
  /
 \multiput {\hbox{$\bullet$}} at "vvbad.dat" 
 \plot  "vvbad.rep" 
 \put {Figure 2} at 4 45
\endpicture
\hfill
}}
\medskip
\quad Nonlinear regression with bad data set
\medskip
{\openup 1pt
\beginexample
 Nonlinear regression with bad data set
 index         value      std dev       1       2       3   
    1   Linf  4.8905e+01 5.9938e+00  1.0000
    2   K     2.1246e-01 1.2076e-01 -0.8923  1.0000
    3   t0   -5.9153e-01 1.4006e+00 -0.6548  0.8707  1.0000
\endexample
}
To invoke the robust regression procedure it is only necessary to
make three changes to the existing code. The template for the
robust regression version of the model can be found in the
file {\tt VONBR.TPL}. 
\beginexample
DATA_SECTION
  init_int nobs;
  init_matrix data(1,nobs,1,2)
  vector age(1,nobs)
  vector size(1,nobs)
PARAMETER_SECTION
  init_number Linf
  init_number K
  init_number t0
  vector pred_size(1,nobs)
  objective_function_value f
  init_bounded_number a(0.0,0.7,2)
PRELIMINARY_CALCS_SECTION
  // get the data out of the columns of the data matrix 
  age=column(data,1);
  size=column(data,2);
  Linf=1.1*max(size);  // set Linf to 1.1 times the longest observed length
  a=0.7;
PROCEDURE_SECTION
  pred_size=Linf*(1.-exp(-K*(age-t0)));
  f=robust_regression(size,pred_size,a);
\endexample
The robust regression routine  needs a parameter, {\tt a}, which  
estimates the amount of contamination by outliers.  
This parameter is declared in the
\PS.
\beginexample
  init_bounded_number a(0.0,0.7,2)
\endexample
\noindent The value of {\tt a}
should be restricted to lie between 0.0 and 0.7 
(See the discussion on robust regression in the AUTODIF user's
manual if you want to know where these numbers come from)
so it is defined to be of type {\tt init\_bounded\_number}.
In general it is not possible to estimate the parameter {\tt a}
until we have done a preliminary fit of the model. To do
this we carry out the the minimization in two phases. During the
first phase {\tt a} should be held constant.  The {\tt 2} in the
declaration for {\tt a} causes {\tt a} to be constant until the second
phase of the minimization.  The second change to the model
involves the default initial value {\tt a}. The default value for
a bounded number is the average of the upper and lower
bounds. For {\tt a} this would be
$0.35$ which is too small. We want to use the upper bound of $0.7$.
This is done by adding the line
\beginexample
  a=0.7;
\endexample
\noindent in the {\tt PRELIMINARY\_CALCS\_SECTION}. Finally we modify
the statement in the {\tt PROCEDURE\_SECTION} 
including the {\tt regression} function to 
\beginexample
  f=robust_regression(size,pred_size,a);
\endexample
\noindent to invoke the robust regression function. That's all there is to
it! These three changes will modify any AD Model builder template
from a nonlinear regression model to a robust nonlinear regression model. 
\vbox{
\medskip
\quad\hbox{
\beginpicture
  \setcoordinatesystem units <.18in,.04in>
  \setplotarea x from 0 to 16.5, y from 0 to 50 
  \axis left label {Size} ticks
    numbered from 0 to 50 by 10 
  /
  \axis bottom label {Age} ticks
    numbered from 0 to 16 by 2 
  /
 \multiput {\hbox{$\bullet$}} at "vonbg.dat" 
 \put {Figure 3} at 4 45
 \plot  "vvrgood.rep" 
\endpicture
\hfill
}}
\medskip
 Robust Nonlinear regression with good data set
\medskip
{\openup 1pt
\beginexample
 index         value      std dev       1       2       3       4   
    1   Linf  5.5707e+01 1.9178e+00  1.0000
    2   K     1.7896e-01 1.9697e-02 -0.9148  1.0000
    3   t0    2.1490e-01 2.0931e-01 -0.5604  0.7680  1.0000
    4   a     7.0000e-01 3.2246e-05 -0.0001  0.0000 -0.0000  1.0000
\endexample
}
The results for the robust regression fit to the bad data set are shown in
figure~4. The estimate for $L_\infty$ is $56.18$ with a standard
deviation of $3.68$ to give a 95\% confidence interval of about
$56.18\pm 7.36$. As expected the both the
parameter estimate and the confidence limit is much less affected 
by the outliers for the robust regression  estimates. The
parameter {\tt a} is estimated to be equal to $0.36$ which indicates
that the robust procedure has detected some rather large outliers.

The results for the robust regression fit to the good data set are shown in
figure~3. The estimates are almost identical to the least-square 
estimates for the same data. This is a property of the robust estimates.
They do almost as well as the least-square estimates when the
assumption of normally distributed errors in the statistical model 
is satisfied exactly, and they do much better than least square estimates
in the presence of moderate or large outliers. You can lose only a little
and you stand to gain a lot by using these estimators. 

\vbox{
\medskip
\quad\hbox{
\beginpicture
  \setcoordinatesystem units <.18in,.04in>
  \setplotarea x from 0 to 16.5, y from 0 to 50 
  \axis left label {Size} ticks
    numbered from 0 to 50 by 10 
  /
  \axis bottom label {Age} ticks
    numbered from 0 to 16 by 2 
  /
 \multiput {\hbox{$\bullet$}} at "vvbad.dat" 
 \put {Figure 4} at 4 45
 \plot  "vvrbad.rep" 
\endpicture
\hfill
}}
\medskip
\quad Robust Nonlinear regression with bad data set
\medskip
{\openup 1pt
\beginexample
 index         value      std dev       1       2       3       4   
    1   Linf  5.6184e+01 3.6796e+00  1.0000
    2   K     1.6818e-01 3.4527e-02 -0.9173  1.0000
    3   t0    6.5129e-04 4.5620e-01 -0.5483  0.7724  1.0000
    4   a     3.6144e-01 1.0721e-01 -0.1946  0.0367 -0.2095  1.0000
\endexample
}
\end
